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ABSTRACT 



O ! In this work we perform a global analysis of the radial velocity curve of the 

^ ! HD 12661 system. Orbital fits that are obtained by the genetic and gradient 

! algorithms of minimization reveal the proximity of the system to the 6:1 mean 

^ ! motion resonance. The orbits are locked in the secular resonance with apsidal 

^ '. axes librating about 180° with the full amplitude ~ (90°, 180°). Our solution 

^ ! incorporates the mutual interaction between the companions. The stability anal- 

^ i ysis with the MEGNO fast indicator shows that the system is located in an 

^ ! extended stable zone of quasi-periodic motions. These results are different from 

O those obtained on the basis of the orbital fit published by Fischer et al. (2003). 

m ■ 
(N ; 

O ! Subject headings: celestial mechanics, stellar dynamics — methods: numerical, 

Q ! N-body simulations — planetary systems — stars: individual (HD 12661) 

P . 1. Introduction 

+3 . 

I The recently discovered planetary systems, HD 12661 and HD 38529 (Fischer et al. 

I 2003), are among about a dozen of multi-planetary systems. Their dynamics are the sub- 
^ I ject of extensive work carried out by many researchers. The knowledge of these systems 

I changes rapidly, as more and more observational data are gathered. Remarkably, most of 
the new multi-planetary systems have a resonant nature. The mean motion resonances 
(MMR) and/or the secular apsidal resonances (SAR) are most likely present in v Andr, 
HD 82943, Gliese 876, 47 UMa, 55 Cnc, and HD 12661 (Fischer et al. 2003, Table 8). 

The orbital fit to the radial velocity (RV) observations of the HD 12661 system has 
been determined by Fischer et al. (2003). They discovered two giant planets, b and c, with 
masses mb — 2.3Mj, and rric — 1.57Mj, revolving around the parent star in elongated orbits. 
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eb ^ 0.35, Cc ~ 0.2 and semi-major axes ab ^ 0.83 au and Qc ^ 2.56 au, respectively. 
These data permit to classify the HD 12661 system as a planetary hierarchical triple system. 
Recently, using the initial conditions provided by Fischer ct al. (2003), Gozdziewski (2003) 
and Lee & Peale (2003), investigated the dynamics of the HD 12661 system. Lee & Peale 
(2003) studied coplanar dynamics in the framework of an analytical secular theory. The 
work of Gozdziewski (2003), mostly numerical, is devoted to the stability analysis of the 
system, in a wide neighborhood of the initial condition, also for non-coplanar systems. The 
results of these papers are in accord. They reveal that the HD 12661 system is close to 
the 11:2 MMR, accompanied by the SAR, with the critical angle ujc — vob librating about 
180°. This type of SAR is the first such case found among the known planetary systems 
(Lee & Peale 2003). The octupole-lcvcl secular theory of these authors explains that the 
planetary system is located in an extended resonance island in the phase space, and it helps 
to understand why the resonance persists in wide ranges of orbital parameters. This is also 
true for non-coplanar systems and wide ranges of planetary masses (Gozdziewski 2003). 

Our dynamical study of the HD 12661 system (Gozdziewski 2003) have been based on 
the first 2-Keplerian fit published by the California and Carnegie Planet Search Team on their 
WEB site^. Recently, an updated fit to the RV data appeared (Fischer et al. 2003, Table 2). 
A preliminary dynamical analysis of this initial condition docs not bring qualitative changes 
to the system dynamics. However, the results of remarkable papers by Laughlin & Chambers 
(2001) and Lee & Peale (2003) inspired us to analyze the RV data again, in the framework 
of the full, nonlinear A'"-body dynamics. Although the RV observations are most commonly 
modeled by additive Keplerian signals, the mutual interaction between giant planets can 
introduce substantial changes to such simple model. A further nuance fiows from the fact 
that the 2-Keplerian fits are mostly interpreted as osculating, astrocentric elements. In fact, 
the fits should be interpreted as Keplerian elements related to the Jacobi coordinates (Lee 
& Peale 2003). These refinements seem to change qualitatively the view of the dynamical 
state of the HD 12661 system, as it can be derived from the current RV data. 



2. Orbital fits 

We used the updated RV observations of the HD 12661 system (Fischer et al. 2003, 
Table 4), containing 86 data points, collected in Keck and Lick observatories. The reported 
uncertainties are ~ 3-6 ms~^ for the Keck data and ^ 7-17 ms^^ for the Lick observations. 

In the first attempt to find the initial condition to the RV data we applied the genetic 
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algorithm (GA). To the best of our knowledge, this method of minimization was used for 

analyzing the RV oft; Andr by Butler et al. (1999) and Stcpinski et al. (2000), the Gliese 876 
system by Laughlin & Chambers (2001), and 55 Cnc by Marcy et al. (2002). Basically, 
the genetic scheme makes it possible to find the global minimum of the function. Our 
experiments confirm remarks of Stepinski et al. (2000) that the method is inefficient for 
finding very accurate best-fit solutions, but it provides good starting points to fast and 
precise gradient methods, hke the Levenberg-Marquardt (LM) scheme (Press et al. 1992). 
These two methods of minimization complement each other. We used publicly available code 
PIKAIA V. 1.2^ by Charbonneau (1995), a great advocate of the genetic algorithms. 

In the simplest case, the minimized function is given by the following formula 



where is the number of observations of the radial velocity V{ti) at moments of time ti, with 
uncertainties a{ti), modeled by a sum of Keplerian signals Vk(^j, p) = Vh{ti, pb) + Vc{ti, p^,), 
and the velocity offset Vq. The fit parameters, p = (Pb,Pc); where Pp = (/Tp, rip, Cp, cUp, Tp), 
for every planet p = b, c, are: the amplitude K^, the mean motion rip, the eccentricity Cp, 
the argument of periastron cUp, and the time of periastron passage Tp. The orbital elements 
emerging from these parameters should be related to the Jacobi coordinates (Lee & Pcale 
2003). Recently, we drew a similar conclusion independently (Gozdziewski et al. 2003). In 
our fit model, the planetary system's reference frame is chosen in such way that the z- 
axis is directed from the observer to the system's barycenter and V{t) is the z-coordinate 
of star velocity. The relation between Kp and the mass of a planet nip, the mass of the 
host star m* and the geometric elements, related to the Jacobi coordinates, is given by 
Kp = (jpapTZpSini, where, for the three-body system, a^, — Gmi,/{m^, + mi,) for the first 
planet, (7c = Gmd (m* -|- mb + mc) for the second companion. For efficiency reasons, the 
code for the dynamical analysis, incorporating the MEGNO indicator, works internally in 
astrocentric coordinates, thus, when it is required, the Jacobi coordinates are transformed 
to these coordinates. In this sense, the astrocentric elements are considered only as a formal 
representation of the initial condition^. Together with , we use a/xF ~ a/x^/^j where 
V = N — Np — lis the number of degrees of freedom and Np is the total number of parameters. 

The PIKAIA code is controlled by a set of 11 parameters. We set the most relevant as 
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foUows: the population number to 256, the number of generations to a rather high value 8000- 
12000, and the variable mutation rate in between 0.005 to 0.07 with crossover probability 
equal to 0.95. The PIKAIA runs, restarted many times, resulted in a number of qualitatively 
similar fits, with a/xE^ — 1-44 and RMS ~ 8.2 ms^^. The RV data have been modeled by 
2-Keplerian signal and one RV offset. We found that, surprisingly, all these solutions, not 
only have a slightly smaller \f)^ and RMS than reported by Fischer et al. (2003) for the 
2-Keplerian model (\/x^ ~ 1.46 and RMS ~ 8.8 ms~^), but also differ from the former fit 
substantially. The periods ratio is very close to 6:1, and not to 11:2! Also, the initial ec ~ 0.1 
is smaller than previous estimates. 

Actually, because the Doppler measurements have been performed in two different ob- 
servatories, we should account for two independent velocity offsets Vq and Vi, different for the 
the Lick and Keck observations. Data that are inhomogeneous in this sense were analyzed 
by Stepinski et al. (2000) and Rivera & Lissauer (2001). Indeed, the GA incorporating the 
RV model with two velocity offsets, makes it possible to obtain a significantly better solution 
(V5d - 1-37 and RMS ~ 7.9 ms'^). 




We verified this fit by determining the quasi-global minimum of as a function of 
(Pc, Cc). We searched for the best fit solutions with these two parameters fixed. Because Pb 
is determined very well, its value, as well as Cb, approximated by the best GA solution for the 
2-Keplerian model, was selected for the starting point in the LM method. Next, varying the 
initial phases (the arguments of pcriastron and the mean anomalies) of the two companions, 
with the step ~ 30°, we calculated the best fit solution with the LM algorithm. The results 
of this experiment are given in Fig. 1. This scan reveals that the best fit solution is localized 
in a fiat minimum. Its parameters agree very well with the GA solution. This test confirmed 
that the GA fit is really global. Finally, we refined this fit by the LM algorithm, driven by 
the full model of the dynamics, incorporating the mutual gravitational interaction between 
the planets. The best GA fit was used as a starting point to the gradient search. The result 
is shown in Table 1, and named the LMl fit. The synthetic RV curve is shown in Fig. 2. 

To estimate the errors of this solution we used, at the first attempt, a method relying 
on a determination of the confidence intervals (Press et al. 1992). We were warned by 
the referees, that this method has drawbacks and the most relevant is that using it we do 
not account for the stellar "jitter". Even for a chromospherically inactive HD 12661 star, 
the uncertainty, (Tjittcr, contributed by the jitter to the RV measurements, can be relatively 
high. The referee, Gregory Laughlin, suggested to us using the estimate (Xjittcr — 5 ms~^, 
based on the data by Saar et al. (1998). Actually, the observations have to be weighted 



by (J = y o"obs + "^fitter' whcre dobs are "pure" instrumental errors. Such joint uncertainty, 
accounted for calculating -^xf, leads to a value ~ 1, thus giving a statistical indication of 
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an adequate model of the RV data. To estimate the parameters errors we synthesized about 
250 sets of "observations". To every original RV measurement we added a Gaussian noise 
with the mean dispersion ctk = 3.4 ms~^, Ul = 7.6 ms~^ for observations gathered by Keck 
and Lick spectrometers, respectively, and the Gaussian noise of the stellar jitter, with the 
dispersion ajitter = 5 ms~^. Next, the Newtonian model of dynamics was fitted to every 
such synthetic data set with the LM algorithm which started from the LMl solution. The 
mean values of the fit parameters are given in Table 1, and called the LM2 fit. Finally, the 
dispersions of these orbital parameters are adopted as the mean uncertainties of both the 
LMl and LM2 fits. In fact, both solutions are the same, with respect to the error bounds. 



3. Stability analysis 

To investigate the dynamical stability of the LMl and LM2 fits we used the fast indicator 
called MEGNO. This technique was advocated by us in a series of recent papers, see, e.g, 
(Gozdziewski et al. 2001; Gozdziewski 2003). The MEGNO indicator is closely related to 
the maximal Lyapunov exponent, but it permits to determine rapidly whether an initial 
condition leads to a regular, quasi-periodic or a chaotic, irregular solution. It is a very 
efficient tool that helps to detect orbital resonances, their structure, and unstable regions in 
the phase space. Looking at a neighborhood of the analyzed initial conditions is a profitable 
way of resolving the question whether the system dynamics are robust to the fit errors. 

The MEGNO tests reveal that both the LMl and LM2 fits are related to quasi-periodic 
motions of the HD12661 system. The MEGNO signature for the LMl fit is shown in Fig. 4c 
and a perfect convergence of the temporal value of MEGNO to 2 indicates a quasi-periodic 
evolution of the planetary system. The MEGNO scan in the neighborhood of this fit, in 
the (oc, ec)-plane, is shown in the left panel of Fig. 3. Actually, both fits are located in 
a relatively extended stable zone, in a proximity of the the 6:1 MMR. The 6:1 MMR is 
separated about 0.2 au from the 11:2 MMR, which has been considered, up to now, as the 
closest MMR neighboring the system in the phase space. The evolution of orbital elements 
in the LMl fit is shown in Fig. 4a,b. 

The stable zone of the 6:1 MMR is relatively very narrow in the range of small tc- To 
illustrate its effect on the motion of the HD 12661 system we changed Oc in the LM2 fit, 
from the nominal value ~ 2.78 au to 2.745 au, still keeping it in the error bound. The 
MEGNO scan, in the (cc, ec)-plane, for such modified initial condition is shown in the right 
panel of Fig. 3, and the Jacobi elements are shown in Fig. 4d,e,f. Note that in both cases 
the SAR with the apsidal lines antialigned in the exact resonance is present, so the system 
still remains in the large libration island found by Lee & Peale (2003). 
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The scan in the (i, AQ)-plane, where AQ — Q^ — is shown in Fig. 5, and it makes is 
possible to estimate the bounds on masses in the systems with the same i, in terms of regular 
and chaotic motions. For coplanar systems, the stable zone extends up to i ~ 25°, but also 
two islands about i ~ 10°, 18° exist. Whether the system can survive during evolutionary 
time scales, in the chaotic regions, can be verified, in general, only by long-term integrations. 

Conclusions 

The study of the RV observations of the HD 12661 star revealed fits which lead to 
significantly different dynamics when compared to the original solution published by Fischer 
et al. (2003). Combined genetic and gradient methods of optimization helped us to find a 
global minimum of ioi the RV model that incorporates the mutual interactions between 
planets. A dynamical analysis shows that the new fits are related to a system close to the 
6:1 MMR, and locked into a deep secular resonance with apsidal lines antiahgned. The 
fits preclude the proximity of the 11:2 MMR, which has been suggested by the previously 
determined initial conditions. In the phase space, the system lies in a relatively extended 
zone of stable, quasi-periodic motions. The closeness of the HD 12661 system to the low- 
order 6:1 MMR adds a new inquiring case to the 3:1 MMR in the 55 Cnc system, the 7:3 
MMR in the 47 UMa system, and the 5:2 MMR of Jupiter and Saturn. It is difficult to 
claim that the HD 12661 system is locked exactly in the 6:1 MMR, but our results suggest 
that very likely it is close to this resonance. Hopefully, the observational window of the 
HD 12661 will cover soon two orbital periods of the outer companion, and the updated set 
of measurements will help to refine the fits again. 

The analysis of the RV data benefits from a global approach, similarly to the stabil- 
ity studies. For planetary systems with large, and not very distant planets, the effects of 
mutual interactions and a proper dynamical interpretation of the orbital fits are vital for 
understanding the dynamics and for resolving all the information hidden in the observations. 
Their analysis requires much care, because omitting even subtle factors can lead to quite a 
different interpretation of the same data set. 
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Table 1: Osculating, astrocentric Keplerian elements of the HD 12661 planetary system for 
the epoch of the first observation (JD=2450831.608). Mass of the host star is 1.07 Mq. 



Fit LMl (A/p=10+2) Mean (LM2) (^^=10+2) 

Parameter Planet b Planet c Planet b Planet c 

mpisini [Mj] ... 2.33 (0.04) 1.69 (0.11) 2.32 (0.04) 1.63 (0.11) 

a[au] 0.822 (0.001) 2.804 (0.11) 0.823 (0.001) 2.781 (0.11) 

e 0.349 (0.016) 0.084 (0.054) 0.343 (0.016) 0.128 (0.054) 

a; [deg] 115.2 (3.9) 292.8 (34.1) 113.8 (3.9) 303.1 (34.1) 

M [deg] 129.4 (4.7) 354.3 (42.9) 132.1 (4.7) 342.5 (42.9) 

Vo (Lick) [ms-i] -7.3 (1.9) -6.3 (1.9) 

Vi (Keck) [ms-^] -12.7 (1.9) -12.2 (1.9) 

1.37 0.93 

RMS [ms-i] .... 7.88 7.44 
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1300 1400 1500 1600 Pc [d] 

Fig. 1. — The minimum values of ^ function of {Pc,ec) for the 2-Keplerian, Jacobi 

model of the RV observations {Np = 10 + 2). The global minimum is marked by the 
intersection of two lines, corresponding to Pc ~ 1630 d, Cc ^ 0.079. See the text for an 
explanation. 
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Fig. 2. — The synthetic radial velocity curve defined by the LMl fit (see Table 1). Time is 
given in JD. Filled circles mark the observations published by Fischer et al. (2003). 
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Fig. 3. — The MEGNO maps in the (Pc, ec)-plane for the initial elements found in this paper. 
The left panel is for the LMl fit, the right panel is for the modified LM2 fit (as explained 
in the text). Zones centered about ^ 2 correspond to regular, quasi-periodic motions of the 
HD 12661 system. The intersection of two lines marks the fitted parameters, the MMRs 
relevant to our discussion are labeled. The resolution of the scan is 120 x 100 data points. 
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Fig. 4. — Tlie evolution of Keplerian orbital elements related to Jacobi coordinates for the 
initial conditions found in this paper. The left column is for the LlVll fit, the right column 
is for the modified L1VI2 fit. Note the presence of the secular apsidal resonance — its critical 
argument — vOc librates about 180° (panels b,e). The modified LM2 fit leads to the 
6:1 IVIIVIR (panel f), while for the LlVll fit the same critical argument circulates. Note the 
stabilizing effect of the 6:1 MMR on the eccentricities of the planets (panels a,d). Panel c is 
for MEGNO signature of the LMl fit and it indicates a quasi-periodic, stable motion of the 
planetary system. Time is given in 10^ yr. 
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Fig. 5. — Stability regions in the plane of the system inclination i vs the relative inclination 
of orbits (expressed through the difference of the longitudes of nodes, AQ = Qc — ^h), derived 
for the LM2 fit. The resolution of the scan is 85 x 90 data points. 



